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We introduce a novel strategy for cosmological Boltzmann codes leading to an increase in speed by 
a factor of ~ 30 for small scale Fourier modes. We (re-)investigate the tight coupling approximation 
and obtain analytic formulae including the octupoles of photon intensity and polarization. Numeri¬ 
cally, these results reach optimal precision. Damping rapid oscillations of small scale modes at later 
times, we simplify the integration of cosmological perturbations. We obtain analytic expressions 
for the photon density contrast and velocity as well as an estimate of the quadrupole from after 
last scattering until today. These analytic formulae hold well during re-ionization and are in fact 
negligible for realistic cosmological scenarios. However, they do extend the validity of our approach 
to models with very large optical depth to the last scattering surface. 


I. INTRODUCTION 

Standard codes such as CMBFAST 0,13 , CAMB iiS 
or CMBEASY lEOli compute the evolution of small 
perturbations in a Friedman-Robertson-Walker Universe. 
The output most frequently used are multipole spec¬ 
tra of the Cosmic Microwave Background (CMB) and 
power spectra of massive particles. These predictions are 
compared to precision measurements of the Cosmic Mi¬ 
crowave Background (CMB) Q and Large Scale Struc¬ 
ture (LSS) [l3- Working in Fourier space, the codes 
evolve perturbation equations for single Fourier fc-modes. 
The simulated evolution starts well outside the horizon 
at early times and ends today. For the CMB, relevant 
scales lie in the range k ~ 10“® ... 1 Mpc“^, while those 
for the LSS extend to higher fc ~ 5 ... 1000Mpc“^. 

Currently, the time needed to evolve a single mode 
is roughly proportional to k. As the spectrum is com¬ 
puted in logarithmic fc-steps, the largest few fc-modes 
tend to dominate the resources needed for the entire cal¬ 
culation. We have analyzed the current strategy to in¬ 
tegrate the perturbation equations and singled out two 
bottlenecks. The first one is the so called tight coupling 
regime (or better: the end of tight coupling). The second 
one are rapid oscillations of relativistic quantities for high 
fc-modes. Roughly speaking, in standard CMBFAST and 
CMBEASY, both regimes contribute equally to the com¬ 
putational cost. This is likely not the case for CAMB, as 
it uses a higher order scheme during tight coupling - a 
solution similar^ to the one we will present later on. 

Our strategy therefore consists of two parts. The first 
one is a revised tight coupling treatment. In this, we 
will make a conceptual change, distinguishing between 
tight coupling of the baryon and photon fluid velocities 
on one hand and the validity of an analytic treatment of 
the photon intensity and polarization quadrupole on the 
other. In essence, our solution extends to the octupole. 


^ To our knowledge, there is no published discussion on higher 
order schemes. There is, howeve r u npublished work by Antony 
Lewis and Constantinos Skordis 
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FIG. 1: The CMB multipole spectrum up to 1 = 4000 for a 
standard cosmological model (solid line). The dashed (blue) 
line shows the relative deviation between a standard CMBEASY 
calculation and one where the switch ending tight coupling 
has been pushed to earlier times and hence better precision. 
The geometric average deviation is ~ 0.3%. The dashed- 
dotted (red) line shows the deviation between such a high 
precision CMBEASY calculation and the new algorithm. With 
the average geometric deviation ~ 0.01% roughly 30 times 
smaller, our new algorithm comes close to the optimal result. 


We thus capture the physics during tight coupling better 
than previously achieved. This leads to a considerable 
increase in accuracy reaching the optimal precision for 
this stage of the computation (see Figure QJ. 

The second part of our solution consists of suppress¬ 
ing unwanted oscillations in the multipole components of 
relativistic particles. In essence, it is the line-of-sight Q 
formulation of all modern CMB codes that allows us to 
do this. As we will see, the oscillations we suppress are 
anyhow unphysical as they perpetuate unwanted reflec¬ 
tions due to truncation effects. In any case, the mod- 
ihcations are such that observational quantities like the 
CMB or LSS are not influenced by our choice. These two 
improvements combined lead to considerably shorter in¬ 
tegration times. Typically, the benefit sets in for modes 
k ^ 0.1 Mpc~^ and increases gradually until reaching fac¬ 
tors of ~ 30 for modes ~ 5 Mpc“^ and higher. For some 
speed comparisons, see Tabled 
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kmax / h, 

Imax 

CMBEASY 

CMBEASY 



(new algorithm) 

(sync, gauge) 

10 Mpc~^ 

no CMB 

1.5s 

10s 

100 Mpe"^ 

no CMB 

4s 

93s 

5Mpc“^ 

2000 

5s 

12s 

lOMpc”^ 

4000 

9s 

25s 


TABLE I: Comparison of speed between the new algorithm 
and the standard synchronous gauge implementation. Ex¬ 
ecution times of CMBFAST are comparable to the standard 
synchronous gauge implementation, but can deviate by a fac¬ 
tor of ~ ^ ■ 2 from CMBEASY depending on the task. The 

Hubble parameter for the model used was h = 0.7. 


II. TIGHT COUPLING REVISED 


At early times, the photon and baryon fluids are 
strongly coupled via Thomson scattering. The mean free 
path between collisions of a photon = anecrr is given 
in terms of the number density of free electrons ne, the 
scale factor of the Universe a and Thomson cross section 
<JT- During early times, Hydrogen and Helium are fully 
ionized, hence rif, oc a~^ and Tc ex. a'^. During Helium and 
Hydrogen recombination, this scaling argument does not 
hold (see Figure 0 ). To avoid these periods we resort 
to the correct value of fc computed beforehand instead 
of using fc = for redshifts z < 10'^. The effect of 
assuming that the scaling holds would however be con¬ 
siderably less than 1% on the final CMB spectrum. 

To discuss the tight coupling regime, let us recapitulate 
the evolution equations for baryons and photons. We do 
this in terms of their density perturbation S and bulk ve¬ 
locity V. For photons, we additionally consider the shear 
a-y and higher multipole moments Mi of the intensity as 
well as polarization multipoles Ei. Our variables are re¬ 
lated to the ones of na by substituting v ^ k ^9. In 
longitudinal gauge, baryons evolve according to 

5b = —kvb + 3(^ (1) 

Vb = —-Vb + clkSb + RT~^{vj — Vb) + k^l;, ( 2 ) 


where R = (4/3)p.y/ pb, the speed of sound of the baryons 
is denoted by and (j) and ^ are metric perturbations. 
By definition, R oc a~^ (provided no baryons are con¬ 
verted to other forms of energy) and at the time of inter¬ 
est, cl (X Tb = (X a~^ (for more detail see e.g. 0 )- 



FIG. 2: Relative deviation of Tc from the naive scaling rela¬ 
tion: [fc. — 2(a/a)Tc]/[2(a/a)Tc] (solid line). We also depict 
the prodnet (dashed line) vs. the scale factor a, which 
compares the mean free path to the expansion rate of the 
Universe. In the cosmological model used, matter radiation 
equality is at aequ = 3 x 10~^ and last scattering defined by 
the peak of the visibility function is at ais = 9 x 10“^. The 
deviation around a = 2 x 10““^ is from Helium recombina¬ 
tion and is practically negligible, because the visibility is still 
small during that period. At later times, however the devia¬ 
tion is due to the onset of Hydrogen recombination and takes 
on substantial values before last scattering. 


Photons evolve according to the hierarchy 


5^ = --kvj + 4:(j) 


= k ( -5.y - ]+ kip + Tc ^{vb - Vj), 


( 3 ) 

( 4 ) 


= M 2 = -Tc ^ I ^M2 + ^E- 


10 


10 


'2 3 

+ ( gt'v - 7-^3 


( 5 ) 


Ml = k 


Mi-I- \ M i+i )-Tc ^Mi, ( 6 ) 


21-1 21 + 3 

where the U-type polarization obeys 


E 2 = -k^Ea - T-^ (^E2 + 


^2 i- 

10 10 , 


El = -t-^Ei 


+k 


21 - 1 


Ei-i- 


y/P + 21-3 
21 + 3 


( 7 ) 


Ei+i .( 8 ) 


The overwhelmingly large value of precludes a 
straight forward numerical integration at early times: 
tiny errors in the propagation of Vb and v-y lead to strong 
restoring forces. This severely limits the maximum step 
size of the integrator and hence the speed of integration. 
Ever since Peebles and Yu fi''st calculated the CMB 
fluctuations, one resorts to the so called tight coupling 
approximation. This approximation eliminates all terms 
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of order from the evolution equations assuming^ tight 
coupling at initial times. Our discussion will closely lean 
on that of ,113, taking a slightly different route. In con¬ 
trast to 1 121 . h owever, we will keep all terms in the deriva¬ 
tion. Like |12|. we start by solving 0 for (vb — Vy) and 
write v-y = Vb + {"^7 — Vb) to get Equation (71) of 


{vb - Vj) = Tc 


Vb + (f’7 


Vb) - k 





(9) 

Substituting Equation 0 , for Vb into this Equation 
one gets Equation (72) of 


( 1 +^) 

Tc 


{vb - Vj) 


- Vb + (t’7 - Vb) 


+ k I ctSb — -S^ -f tj. 


( 10 ) 


where we have used = — fc^. We could stop here, 
however it is numerically better conditioned to write 
^ (a) ^here ^Vb is obtained from solv¬ 
ing Equation for ^Vb- This expression for is 

then plugged into Equation m to yield the final result 
for the slip (denoted by V) 


V = {vb- U7) = 


1 + R 


1 + R 


-{vb - Vj) 


1 


- Vb + {Vj -Vb) + k[ -Sj - 2(7^ -I- 'll; 


k ( CgSb — T '^7 T *^7 


, Tr 

1 + 2 -- 


al + R 


(15) 


or alternatively, at times when the scaling of Tc holds. 


Deriving the LHS of this Equation yields 

(1 + i?), 


Ihs = 


-{vb - Vj) 
(1 + R) 


(vb - U7) 

d R 1 + Rfc 


a Tc 


( 11 ) 


('Vb - U 7 ) - ^ ^^ -(Vb - U 7 ), (12) 

Tc a 


where the last line holds provided the assumed scaling 
of Tc is correct (see also Eigure El). All in all, deriving 
Equation m with respect to conformal time yields 


(1 + i?) 


('Vb - 1)7) - 


d R 
a Tc 


= (U7 - Vb) - 


a 


1 + RTc 

Tc Tc 
2 


(vb - Vj) 


1^7 *^0 J 

a 


Cl . 

Vb - Vb 

a 

1 


+ k ( CgSb CgSb — “ 7(^7 H" 


(13) 


Multiplying Equation © by I to substitute ^■Ub in m, 
we get 


(1 -I- i?) , 

- (Vb - Vj) = 


d R 1 + Rfc 
a Tc 


+ (v^ - Vb) - Vb -I- 2 - Vb - 2 -cikSb 


Tc Tc 
2 


(vb - Vj) 


a 

1 , 
4 


-I- k ( c^Sb — -^7 -I- CT 7 ) -I- -kip (14) 


^ There is no restoring force left, as we will see. Any error in 
the approximation is therefore amplified over time. One could, 
in principle retain a fraction of the restoring force to eliminate 
small numerical errors. However, this is not necessary in practice 
and we therefore will not discuss this possibility further. 




- -Vb + (v^ - Vb) + k [ i(57 - 2(7^ -I- ip 


1 + R 


+ k ( Cf,Sb T <77 


1-f 2 


2 •'7 
d Tc 


al + R 


(16) 


This Equation itT^ l or more obviously (EU) is essentially 
Equation (74) of [l2| up to some corrections. Having kept 
all terms, we note that our Equation E3 is exact. To 
obtain Equations of motion for Vb and ^7 during tight 
coupling, we plug our result for (u{, — U 7 ), Equation 1151) 
into the RHS of Equation © and this in turn into the 
RHS of Equations © and 0 . This yields 


Vb = 


^ kc^Sb - -Vb ) -I- kip 


1 + R 


R 


1 + R 


/C ( ^ (7.y 


V 


^ , /^Ir 1 , , 

-K ( -d.y — CTj I + kip 


1 + R V4 


( 17 ) 


Up to now, we have made no approximations. Conceptu¬ 
ally, we would like to separate the question of tight cou¬ 
pling for the velocities v.y and Vb from any approximations 
of the shear which we make below. As far as the tight 
coupling of the velocities and hence the slip V is con¬ 
cerned, our approximation is to drop the term (1).^ — Vb)- 
We reserve the expression ’tight coupling’ for the validity 
of our assumption that ( 1)7 — Vb) can be neglected in the 
slip V. As a criterion, we use kvc < ^ for the photon 
fluid. When this threshold is passed, we use Equation 0 
to evolve the photon velocity. Likewise, for the baryons, 
we use max(fc, ^)tc/R < Again, when this limit is 
exceeded, we switch to Equation ©. In any case, we 
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FIG. 3: The quadrupole M 2 obtained by a full numerical 
evolution for a mode of fc = IMpc”*^ (dotted line). The 
solid (blue) line depicts the deviation of our analytic result, 
Equation from the numerical value. For this mode, we 
normally switch to the full numerical evolution at r = 65 Mpc 
when the analytic estimate still holds very well. 


switch off the approximation At = 30 Mpc before the 
first evaluation of the CMB anisotropy sources (see be¬ 
low). For a A — CDM model, this is at r « 200 Mpc. 

To obtain high accuracy during tight coupling, it is 
crucial to determine Not so much for the slip dia, 
but more so for the Equations of motion 03: the shear 
reflects the power that is drained away from the velocity 
in the multipole expansion. This leads to an additional 
damping for photons. For the shear, we distinguish two 
regimes: an early one, where we use a high-order analytic 
approximation and a later one in which the full multipole 
equations of motion are used. 

Since Tc 1 at early times, one gets from multiply¬ 
ing © by Tc that Mi « (krc) Mi-il/{21 — 1). Hence, 
higher multipoles are suppressed by powers of krc. Ap¬ 
proximating this situation hy M3 = E3 = M4 = E4 = 0 
in Equations and m, we get 

Ms = ^ {kTc)M2 
5 

Es = ■^{kT,)E2. ( 18 ) 

Likewise, we obtain a leading order estimate of the 
quadrupoles by temporarily setting M 2 = E 2 = Q, 

= l{kT,)v, ( 19 ) 

( 20 ) 

Inserting Equations CHI into the quadrupole Equations 
© and 0 and using M 2 = M 2 °' and E 2 = ^ 2 °' an 
estimate for the derivative, we get the desired expression 
for the shear 


-cr^ = M 2 = -kTcV^ 


1 _ ^ {kr^f 
70 ^ ’ 


-^TcM^ 2 ' 
6 ^ 


( 21 ) 


which is precise to order Tc and {kTc)^ (see also Eigure 
©. The inclusion of the octupole reduces the power of 
M 2 as expected. 

In practice, we use A42°’ to calculate the slip V* ° to 
leading order. This in turn is used to calculate From 

vi^°, we get M 2 °' which in turn is needed to obtain the 
accurate value of A42 according to Equation (ETll . The 
difference A 7 VI 2 = M 2 — M 2 °' is then used to promote 
yi.o 1 ) g^g gg ^Lo. Finally, having A42 and 

V at hand, we get hf, from Equation 03- 

When this approximation breaks down (sometimes 
long before tight coupling ends), we switch to the full 
multipole evolution equations. Tight coupling is appli¬ 
cable for kTc <C 1. Equation ED on one hand goes to 
higher order in kTc, namely, as M 2 °' is already of order 
(kTc), our results incorporates quantities up to (krc)^. In 
terms of Tc alone, however. Equation 113 is accurate to 
order Tc (krc) only. Hence, when Tc reaches ~ 10“^ Mpc, 
our analytic expression is not sufficiently accurate any¬ 
more. This signals the breakdown of our assumption that 
M 3 = £3 = 0 (and likewise for higher multipoles). Luck¬ 
ily, it is not critical to evolve the full multipole equations 
even when is still substantial. This is in strong con¬ 
trast to the coupled velocity equations which are far more 
difficult to evolve at times when the analytic quadrupole 
formulae breaks down. In essence, distinguishing be¬ 
tween tight coupling and the treatment of the quadrupole 
evolution is the key to success here. 


III. A CURE FOR RAPID OSCILLATIONS 

While the gain in speed from the method described in 
the last section is impressive, high A:-modes would still 
require long integration times. To see this, one must 
consider the evolution of the photon and neutrino multi¬ 
pole hierarchies.^ Our discussion is aimed at small scale 
modes which are supposed to be well inside the horizon, 
i.e. fcr ^ 1. 

Before last scattering, (fcrc) <C 1 and Mi cx {kTcf~^ 
for I > 1 and so the influence of higher multipoles 
on Sj and Vj may be neglected to first order. In the 
small scale limit that we are interested in, Sj and Vj 
are oscillating according to ^ cos(c2'A:t) and vi, ~ 
sin(c2'fcT). As the speed of sound of the photon-baryon 
fluid is c] « y^l/3, we encounter oscillations with pe¬ 
riod At ft! {2Tr)/{kcj) ft 11/fc. Estimating the time of 
last scattering with tis ft 280 Mpc, we see that a mode 
will perform tis/At ft 25fcMpc oscillations until last 
scattering. Yet, there are many more oscillations after 
last scattering which we turn to now. After last scatter¬ 
ing, is negligible and the multipole hierarchy of pho¬ 
tons effectively turns into recursion relations for spheri- 


® We include the monopole <5-, and dipole Vj here. 
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FIG. 4: The quadrupole as a function of conformal time r 
for a mode of k/h — 0.5Mpc“^ and h = 0.7. The multipole 
expansion for photons and neutrinos has been truncated at 
Imax = 200 (solid line) and Imax = 8 (dashed line) respec¬ 
tively. In the case of Imax = 8, power reflected back from the 
highest multipole Imax renders the further evolution of the 
quadrupole unphysical. Indeed the magnitude of the physi¬ 
cal oscillations are much smaller than the reflected ones. For 
Imax = 200, reflection effects dominate the evolution from 
T ~ 1300 Mpc on. In both cases, the effect of shear of real¬ 
istic particles on the potentials (j) and ip is negligible by the 
time the truncation effects set in. 


cal Bessel functions. The same is true for neutrino multi¬ 
poles which roughly evolve like spherical Bessel functions 
from the start. Spherical Bessel functions have a lead¬ 
ing order behavior similar to jiikr) oc (A:t)“^/^ sin(A:r) 
for fcr ^ 1 and kr > 1. The period is then given by 
At = ( 27 r)/fc. The time passed from last scattering to 
today, is tq — tjs si tq « 14000 Mpc for current cosmo¬ 
logical models. So we encounter ~ tq/At- « fcro/( 27 r) 
=2200 X k Mpc oscillations. Numerically, each oscilla¬ 
tion necessitates ~ 20... 40 evaluations of the full set 
of evolution equations. We therefore estimate a total of 
^ 6 X lO'^ X k Mpc evaluations induced by the oscillatory 
nature of the solution. So a mode k = 5Mpc“^ needs 
~ 3 X 10^ evaluations - a substantial number. 

Since the introduction of the line-of-sight algorithm, 
what one really needs for the CMB and LSS are the low 
multipoles up to the quadrupoles. In fact, the sources for 
temperature and polarization anisotropies are given by 


Sj, = 

'1 


<p + ip 


Vb 


^ + -C 

k P 

C 




and 


■ g —Sry —-—!“((('+ '!/')T7:T 


Se = ^-C (/c[to -t]) ^ 


3 


2 k^ 


C 


( 22 ) 


(23) 


Here, g = Kexp(K(T) — k{tq)) is the visibility with 
k = T~^ the differential optical depth and C = (A42 — 
•\/6£’2)/10 contains the quadrupole information. The 
role of higher multipole moments is therefore reduced to 


FIG. 5: Photon density contrast S~f (upper solid [blue] line), 
(lower solid [black] line) and quadrupole M 2 = (5/2)(j^ 
(dashed dotted [indigo] line) as a function of conformal time 
T before and after re-ionization at r « 2800 Mpc. The [green] 
upper dashed line is the analytic estimate for Sry, Equation 
iHZli and the lower [red] dashed line is the analytic estimate for 
Vj, Equation The analytic estimate of <5-, falls almost on 
top of the correct numerical result. Please note the different 
scales for <5-, and and A42 respectively. The quadrupole 
is roughly of the same order as v-y. The mode shown is for 
k/h = 5Mpc“^ where h — 0.7 and the optical depth to the 
last scattering surface is Topt = 0.3. Please note that we 
truncated the multipole hierarchy at sufficiently high Imax = 
2500. With insufficient Imax, rapid unphysical oscillations of 
considerably higher amplitude would be present. 


draining power away from S^, Vj and M 2 and E 2 (and 
likewise for neutrinos). As the oscillations are damped 
and tend to average out, it suffices to truncate the mul¬ 
tipole hierarchy at low I ~ 8... 25 in the line-of-sight 
approach. This is one of the main reasons for its su¬ 
perior speed. Truncating the hierarchy, though leads to 
unwanted reflection of power from the highest multipole 
Imax- As one can see in Figure 0 the power reflected 
back spoils the mono frequency of the oscillations. At 
best, the further high frequency evolution of the multi¬ 
poles is wrong but negligible, because the oscillations are 
small and average out. This is indeed the case in the 
CMBFASt/camb/cmbeasy truncation. 

We will now show that the overwhelming contribution 
from Sj and C (and its derivatives) of some small scale 
mode k > 10“^ Mpc“^ towards CMB fluctuations comes 
from times before re-ionization. To do this, let us find an 
analytic approximation to the photon evolution after de¬ 
coupling and in particular during re-ionization. Without 
re-ionization, and neglecting AI 2 as well as using </> re ip 
and (p re Q, the equation of motions and 0 can be 
cast in the form 


^7 = -^fc^Q<57 + V'), (24) 

which has the particular solution 

S-y = —Alp. 


( 25 ) 
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As the oscillations of and higher multipoles are 
damped roughly oc we see that to good ap¬ 

proximation, (5-y = —4-0 after decoupling (and before re¬ 
ionization) and all higher moments vanish. 

During re-ionization, reaches moderate levels 

again. As Vb has grown substantial during matter dom¬ 
ination, the photon velocity starts to evolve towards 
Vb- Any increase in magnitude of w-,,, is however swiftly 
balanced by a growth of J-y according to Equation Q. So 
roughly speaking, during re-ionization, we may approxi¬ 
mate 






(26) 


where we omit the tiny term and (a bit more wor¬ 

risome) M. 2 - Hence, during re-ionization, the particular 
solution to the equation of motion is 

(27) 

IvTq 

This approximation holds well (see FigureEl) and oscilla¬ 
tions on top of it are again damped and tend to average 
out. Deriving the above m, one gets 


Vy 


k \ kTc kTc Tc) 


(28) 


Please note that during the onset of re-ionization, fc = 
2 -Tc does not hold and it depends on the details of the re¬ 
ionization history to what peak magnitude will reach. 
Both CMBFAST and CMBEASY implement a swift switch 
from neutral to re-ionized and it is likely that both serve 
as upper bounds on any realistic contribution of higher 
k modes towards the CMB anisotropies at late time. In 
other words: as the effects are negligible for the currently 
implemented re-ionization history, they will be even more 
so for the real one. Going back on track, we give an 
estimate for the amplitude of Af 2 : assuming Ad; « 0 and 
~ 0, one gets from the equations of motion 
that neighboring multipoles A4; are of roughly the same 
amplitude. So the amplitude of Ai 2 and hence that of 
the shear is related to that u-y, i.e. we find the bound 


max(|cr-y|) ~ max(|r!.y|), (29) 


where it is understood that the maximum is taken of 
full oscillations. After radiation domination, the metric 
potential ip is given by 


c: Pc^c 
2M|fc2’ 


(30) 


where Mp is the reduced Plank mass, pc is the energy 
density of cold dark matter and 6 c is its relative density 
perturbation. For modes that enter the horizon during 
radiation domination, Sc is roughly independent of scale 
(we omit the overall dependence on the initial power spec¬ 
trum in this argument). Hence, ip cx during matter 



FIG. 6 : Cold dark matter power spectrum using the old gauge 
invariant implentation (dashed line) and the new strategy 
in gauge invariant variables (thin solid line). The density 
contrast shown is the gauge invariant combination = 

giongzt. _ 20 mean deviation between the curves is 

« 0.02%. To guide the eye, we also depict the synchronous 
gauge power spectrum [thin gray dotted line]. The difference 
at large scales is due to gauge ambiguities. Again, we used 
h = 0.7. 


domination and we see that ip ^ Q and so <5.y, —*■ 0 ac¬ 
cording to Equation m- Provided that Tc/tc remains 
reasonable, v,y and hence M .2 and E 2 will remain negli¬ 
gible as well during re-ionization and afterwards. 

For the LSS evolution, neglecting the shear is a good 
approximation because Einstein’s Equation gives 

[p^Al2 -f = Mpfc^((/)- V'), (31) 

where N 2 is the neutrino quadrupole. As oc a~'^, 
the difference of the metric potentials vanishes for small 
scale modes, i.e. at least 

{(p — Ip) (X {ka)~'^, (32) 

where we have neglected the decay of the quadrupoles 
M 2 and A /2 which give an additional suppression (see 
also Figure EJ. 

As the effect of 6 .y and AT 2 and E 2 at late times for 
small scale modes can be neglected (or very well approx¬ 
imated in the case of 6-y), we see that there is really no 
need to propagate relativistic species at later times. The 
key to our final speed up is therefore to avoid integrating 
these oscillations after they have become irrelevant. We 
do this by multiplying the RHS of equations 0-IHl) as 
well as the corresponding multipole evolution equations 
for relativistic neutrinos by a damping factor P. Defin¬ 
ing X = kr, we employ P = {I — tanh([a; — Xc\/w)}/2 
with the cross over Xc = max(1000, fer^iZute), where 
a{Tdiiute) = ^o-equ and ttequ IS the scale factor at matter- 
radiation equality. This later criterion ensures that the 
contribution of relativistic species to the perturbed en¬ 
ergy densities is negligible: from equality on, 6 c oc a, 
whereas 6 -y decays and pdPrei oc a~^ so at least 

6cPc ■ S-^p-^ cx CL 


(33) 
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and similar arguments hold for neutrinos. Hence, from 
Tdiiute on, one can safely ignore this contribution. The 
former criterion Xc < 1000 ensures that oscillations have 
damped away sufficiently. The cross-over width w is 
rather uncritical. We used w = 50 to make the transition 
smooth. Typically, TdUute ~ 400 Mpc and one therefore 
has to follow only a fraction of TdUute/TQ oscillations as 
compared to the standard strategy. This corresponds to 
a gain in efficiency by a factor TQ/rdUute ~ 30. 

To compute the sources St and use the expres¬ 

sions 


5^ 

C 

C 

c 


Y^^numeric. 

i 0^ 


-4(i-r) 




-pQnumeric. 

Y^QUumeric. 

YQnumeric. 


Vb 

kxc 


(34) 

(35) 

(36) 

(37) 


which interpolate between the numerical value before T- 
damping and the analytic approximations, Equation m 
and C = 0. Setting C = 0 is an approximation to the 
small value of the quadrupoles averaged over several os¬ 
cillations. 

For general dark energy models with rest frame speed 
of sound Cg > 0 of the dark energy fluid, the dark energy 
perturbations well inside the horizon oscillate with high 
frequency. In this case, one needs to suppress the damped 
oscillations of the dark energy fluid perturbations much 
like those of photons to achieve faster integration. 


IV. CONCLUSIONS 

We have improved the integration strategy of modern 
cosmological Boltzmann codes. As a first step, we made 


a conceptual distinction between tight coupling of the ve¬ 
locities Vj and Vb and the validity of analytic estimates 
for the intensity and polarization quadrupole. Doing so 
allowed us to switch to the full numerical evolution later. 
The inclusion of shear at early times lead to an increase 
in precision. In the second part of our work, we inves¬ 
tigated the behavior of photons after decoupling. We 
found analytic approximations for both 6 ^ and as well 
as a bound on the shear a-y. The contributions of photons 
and neutrinos towards CMB anisotropies can be well ap¬ 
proximated by using these analytic estimates of Sy and 
ay for small scale modes deep inside the horizon. In fact, 
for an optical depth Topt ^0.2, late time effects of pho¬ 
tons on the CMB anisotropy sources St and Se may be 
neglected altogether on small scales. We introduced a 
smooth damping of high frequency oscillations of photon 
and neutrino multipoles. The damping effectively freezes 
their evolution well inside the horizon. All in all, our 
strategy leads to a gain in efficiency of up to factor ^ 30 
and comes close to optimal accuracy for both the CMB 
and LSS. 
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